function map = glm_SingleStudy_tMap(hfile, c, mapopts)
% GLM::SingleStudy_tMap  - calculate a t contrast map
%
% FORMAT:       map = glm.SingleStudy_tMap([c, mapopts])
%
% Input fields:
%
%       c           NxC contrast vector (default: full model and main eff)
%       mapopts     structure with optional fields
%        .interp    mesh-based interpolation (default: true)
%        .srf       surface file, required for interpolation
%
% Output fields:
%
%       map         MAP/VMP/SMP object with C maps

% Version:  v0.7c
% Build:    7120508
% Date:     Dec-05 2007, 08:21 AM CET
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'glm')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if bc.ProjectTypeRFX > 0 || ...
    bc.NrOfStudies ~= 1
end
numtp = bc.Study(1).NrOfTimePoints;
numpred = bc.NrOfPredictors;
nval = numtp - numpred;
szmap = size(bc.GLMData.MCorrSS);
if nargin < 2 || ...
   ~isa(c, 'double') || ...
    isempty(c) || ...
    size(c, 1) > numpred || ...
    any(isinf(c(:)) | isnan(c(:)))
    c = [[eye(numpred - 1);zeros(1, numpred - 1)], [ones(numpred - 1, 1);0]];
elseif numel(c) == size(c, 2)
    if numel(c) > numpred
        c((numpred + 1):end) = [];
    end
    c = real(c');
else
    c = real(c);
end
nummaps = size(c, 2);
if size(c, 1) < numpred
    c = [c; zeros(numpred - size(c, 1), size(c, 2))];
end
if nargin < 3 || ...
   ~isstruct(mapopts) || ...
    numel(mapopts) ~= 1
    mapopts = struct;
end
ipo = true;
if isfield(mapopts, 'interp') && ...
   ~isempty(mapopts.interp) && ...
   (isnumeric(mapopts.interp) || islogical(mapopts.interp))
    ipo = ipo && mapopts.interp(1);
end

% generate map according to type
iXX = bc.iXX;
switch (bc.ProjectType)
    
    % currently unsupported types
    case {0, 2}
        error( ...
            'BVQXfile:NotYetSupported', ...
            'FMR/MAP and SRF/SMP generation not yet supported.' ...
        );
    
    % VTCs
    case {1}
        
        % create VMP
        bbox = aft_BoundingBox(hfile);
        map = newnatresvmp(bbox.BBox, bc.Resolution, ones(1, size(c, 2)));
        mapc = bvqxfile_getcont(map.L);
        
        % calculation of SE map
        semap = sqrt( ...
            (1 - (bc.GLMData.MultipleRegressionR .^ 2)) .* ...
            bc.GLMData.MCorrSS / nval);
        semap(semap == 0) = Inf;
        
        % generate tmaps
        for cc = 1:size(c, 2)
            
            % map options
            mapc.Map(cc).DF1 = nval;
            mapc.Map(cc).DF2 = 0;
            mapc.Map(cc).BonferroniValue = bc.NrOfVoxelsForBonfCorrection;
            mapc.Map(cc).UseRGBColor = 0;
            
            % calculate t stats
            mapc.Map(cc).VMPData = reshape(glmtstat(c(:, cc)', ...
                bc.GLMData.BetaMaps, iXX, semap), szmap);
            
            % set name
            mapc.Map(cc).Name = ...
                sprintf('Contrast:%s', sprintf(' %d', c(:, cc)'));
        end
        
     % unknown types
    otherwise
        error( ...
            'BVQXfile:UnknownOption', ...
            'Unknown GLM ProjectType: %d', ...
            bc.ProjectType ...
        );
end

% set back to memory
bvqxfile_setcont(map.L, mapc);
